NSite = 4;
Source = [ 1 N/4; 1 N];
Target = [N/2+1 3*N/4; 1 N];
%Source = [ 1 N/4; 1 N/4];
%Target = [N/2+1 3*N/4; N/2+1 N*3/4];
Source = [ 1 8; 1 8];
Target = [1 8; 9 16];
% Source = [ 1 N/4; ...
%	   1 N/4];
% Target = [N/2+1 3*N/4; ...
%          3*N/4+1 N];
beta = 500;

% display('Case 1');
% mu = SD(NSite);
% RhoExact = real(2d0 * EF * diag( 1./ (exp(beta*(D-mu))+1d0) ) * EF');
% S = testrank( RhoExact, Source, Target );
% semilogy(S, 'b');
% hold on

% clear RhoExact

% display('Case 2');

% mu = SD(2^2*NSite);
% RhoExact = real(2d0 * EF * diag( 1./ (exp(beta*(D-mu))+1d0) ) * EF');
% S = testrank( RhoExact, Source, Target );
% semilogy(S, 'r');

% clear RhoExact

display('Case 3');

mu = SD(2^4*NSite);
RhoExact = real(2d0 * EF * diag( 1./ (exp(beta*(D-mu))+1d0) ) * EF');
S = testrank( RhoExact, Source, Target );
semilogy(S, 'k');

clear RhoExact


Green = inv(beta * (H - mu*eye(N*N)) + 1i*pi*eye(N*N));
S = testrank( Green, Source, Target );
hold on
semilogy(S, 'r');

